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Abstract 

It is well-known that in fluid dynamics an alternative to customary direct solution methods 
(based on the discretization of the fluid fields) is provided by so-called particle simulation methods. 
Particle simulation methods rely typically on appropriate kinetic models for the fluid equations 
which permit the evaluation of the fluid fields in terms of suitable expectation values (or momenta) 
of the kinetic distribution function /(r, v, t), being respectively r and v the position an velocity of a 
test particle with probability density /(r,v,i). These kinetic models can be continuous or discrete 
in phase space, yielding respectively continuous or discrete kinetic models for the fluids. However, 
also particle simulation methods may be biased by an undesirable computational complexity. In 
particular, a fundamental issue is to estimate the algorithmic complexity of numerical simulations 
based on traditional LBM's (Lattice-Boltzmann methods; for review see Succi, 2001 [l]). These 
methods, based on a discrete kinetic approach, represent currently an interesting alternative to 
direct solution methods. Here we intend to prove that for incompressible fluids fluids LBM's may 
present a high complexity. The goal of the investigation is to present a detailed account of the 
origin of the various complexity sources appearing in customary LBM's. The result is relevant to 
establish possible strategies for improving the numerical efficiency of existing numerical methods. 

PACS numbers: 47.10.ad,05.20.Dd 
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I. INTRODUCTION: BASIC MOTIVATIONS 



Fluid Dynamics represents, is a sense, one of the unsolved problems of contemporary sci- 
ence. In fact, although fluid equations, i.e., the complete set of differential equations which 
are associated to a fluid are in many cases well known (and therefore to be considered as phe- 
nomenological equations), a complete knowledge of their solutions is not achievable. This 
justifies the constant efforts placed in Computational Fluid Dynamics (CFD) to develop 
more efficient numerical simulation methods, particularly suitable for software implementa- 
tion on parallel supercomputers and able to simulate the small-scale dynamics of complex 
fluid systems. This explains the increasing role of CFD both for theoretical research and 
applied sciences, in particular for the development of industrial applications. CFD has the 
goal of determining the numerical solutions of suitable fluid equations, associated to a pre- 
scribed fluid, which are expressed in terms of appropriate physical observables, the so-called 
fluid fields. For neutral isothermal fluids the fluid equations are the incompressible Navier- 
Stokes equations (INSE), which are represented respectively by the Navier-Stokes equation, 
advancing in time the fluid velocity V(r, t), the Poisson equation which, at a given time, 
determines the fluid pressure p(r, t) and the isochoricity condition which forces the fluid 
velocity to result everywhere divergence-less. 

A basic aspect of numerical algorithms is their computational complexity, namely the 
number of (elementary) discrete logical operations which must be performed in a given time 
interval to carry out a prescribed calculation. In the case of CFD this is measured by the 
number of discrete operations required to advance in time the fluid fields in a prescribed 
finite time interval. The computational complexity of numerical solution methods for INSE 
depends obviously both on the "physical" properties of the fluid, i.e., on the characteristic 
spatial and time scales of the fluid fields (L, T), as well as on the the accuracy with which the 
solutions are actually determined. In particular, the numerical methods used in CFD and 
CMFD for incompressible fluids (see for example J. Kim, 1999 j^j]), exhibit typically a high 
computational (or algorithmic) complexity, although due to different reasons. Among the 
most popular are the so-called direct solution methods, which are based on 
of the fluid fields (see, for example J. Kim and P. Moin, 1979 and 1985 3|, 
(ll; J. Ferzinger and M. Peric, 1996 [6]; S. Turek, 1999 7]). Especially for LES (large eddy 
simulations) of fluid problems with complex boundaries and in the presence of strong fluid 
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turbulence these numerical methods are difficult to implement on supercomputers based 
on parallel architectures, mostly due to the difficulty of solving the fluid equation which 
advances on time the fluid pressure (i.e., the Poisson equation). In the worst cases (i.e., 
when the turbulence scale length becomes comparable to the size of the spatial cells and 
is sufficiently widespread in the domain of the fluid), the latter exhibits a computational 
complexity proportional to A 7 , where the exponent 7 is typically larger than 1 and in the 
worst case can be 7 = 2 and N is the number of nodes (or sets) in which the fluid is 
discretized (and which measures the "size" the numerical simulation). This phenomenon 
may result, in actual numerical experiments on incompressible fluids, in the appearance of 
a computational bottleneck which effectively limits the size of numerical simulations. Since 
all fluids with sufficiently small Mach number behave as incompressible, it is obvious that 
this phenomenon represents a major issue in computational fluid dynamics and a formidable 
obstacle to the development of large scale numerical simulations relevant for actual industrial 
applications. 

A key aspect of CFD lies also in the search of efficient parallel numerical algorithms, 
namely algorithms which can be represented in terms of independent or weakly related 
sub-algorithm to be handled independently by a prescribed set of processors working in 
parallel. It is well known that customary, most algorithmically-efficient, direct-solution 
methods adopted in CFD, become highly inefficient for parallel processing. This is due 
to the Poisson equation, an elliptic PDE which results difficult to solve numerically (the 
best numerical methods turn out to be the least efficient for parallel processing). This has 
motivated the search of numerical algorithms which are able to determine the fluid pressure 
without actually solving numerically the Poisson equation. 



II. LATTICE-BOLTZMANN METHODS 



The investigation of lattice Boltzmann methods (LBM) has drawn considerable attention 
in the last few years (for a review see 8j). The simplicity of the algorithms, along with 
the ability to determine the fluid pressure without actually requiring the explicit numerical 
solution of the Poisson equation, have allowed LBM's to emerge as possible alternatives to 
traditional CFD approaches based on the direct discretization of the relevant fluid equa- 
tions (direct solution methods). A specialized interesting area of investigation concerns, 
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in particular, the hardware implementation of LB algorithms for reconfigurable computing 



9|, [lOj, a technique which appears especially relevant for applications to real-time simula- 
tions and process optimization in industrial fluid dynamics. In addition, these approaches 
exhibit a natural parallelism and - at the same time - result highly scalable on an array of 
independent processing units working in parallel. In essence, this is because the elementary 
process (represented by the evolution of the fluid fields) is reduced to the time evolution of 
a suitable set of weakly correlated test particles, each one represented by a discrete kinetic 
distribution function. In customary LBM's typically the fluid is treated as weakly compress- 
ible, with asymptotically small deviations from the condition of true incompressibility. This 
is realized by actually replacing INSE with a modified set of fluid equations. In particular 
the pressure is determined in this case by an equation of state, while the fluid velocity is 
advanced in time by means of the kinetic distribution function. This permits to avoid the 
computational complexity of Poisson equation and to obtain, at the same time, a highly 
parallelizable approximate solution method. For this reason LB approaches are currently 
considered as a promising computational tool for massive fluid-dynamics numerical simula- 
tions in a wide spectrum of applications, ranging from isothermal incompressible fluids, to 
thermal, compressible fluids, to MHD dynamics in fusion plasmas. 

III. DIFFICULTIES WITH LBM'S 

Despite the significant number of theoretical and numerical pap ers appeared in the lit- 



erature in the last few years, the lattice Boltzmann method 
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16| - among 



many others available in CFD - is probably the one for which a complete understanding is 



not yet available. Although originated as an extension of the lattice gas automaton [171. 118] 
or a special discrete form of the Boltzmann equation Jjj] , several aspects regarding the very 
foundation of LB theory still remain to be clarified. Consequently, also the comparisons and 
exact relationship between the various lattice Boltzmann methods (LBM) and other CFD 
methods are made difficult or, at least, not yet well understood. Needless to say, these com- 
parisons are essential to assess the relative value (based on the characteristic computational 
complexity, accuracy and stability) of LBM and other CFD methods. In particular the rela- 
tive performance of the numerical methods depend strongly on the characteristic spatial and 
time discretization scales, i.e., the minimal spatial and time scale lengths required by each 
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numerical method to achieve a prescribed accuracy. On the other hand, most of the existing 
cnowledge of the LBM's properties originates from numerical benchmarks (see for example 



20, 
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22i|). Although these studies have demonstrated the LBM's accuracy in simulating 



fluid flows, few comparisons are available on the relative computational efficiency of the LBM 
and other CFD methods 
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23( | . The main reason [of these difficulties] is probably because 



current LBM's, rather than being exact Navier-Stokes solvers, are at most asymptotic ones 
(asymptotic LBM's), i.e., they depend on one or more infinitesimal parameters and recover 
INSE only in an approximate asymptotic sense. The motivations of this work are related to 
some of the basic features of customary LB theory representing, at the same time, assets and 
weaknesses. One of the main reasons of the popularity of the LB approach lays in its simplic- 
ity and in the fact that it provides an approximate Poisson solver, i.e., it permits to advance 
in time the fluid fields without explicitly solving numerically the Poisson equation for the 
fluid pressure. However customary LB approaches can yield, at most, only asymptotic ap- 
proximations for the fluid fields. This is because of two different reasons. The first one is the 
difficulty in the precise definition of the kinetic boundary conditions in customary LBM's, 
since sufficiently close to the boundary the form of the distribution function prescribed by 
the boundary conditions is not generally consistent with hydrodynamic equations. The sec- 
ond reason is that the kinetic description adopted implies either the introduction of weak 
compressibility Q, [ll|, Q, 14, Q, [3] or temperature (3] effects of the fluid or some sort of 
state equation for the fluid pressure [3] . These assumptions, although physically plausible, 
appear unacceptable from the mathematical viewpoint since they represent a breaking of the 
exact fluid equations. Moreover, in the case of very small fluid viscosity customary LBM's 
may become inefficient as a consequence of the low-order approximations usually adopted 
and the possible presence of numerical instabilities (see below). These accuracy limitations 
at low viscosities can usually be overcome only by imposing severe grid refinements and 
strong reductions of the size of the time step. This has the inevitable consequence of raising 
significantly the level of computational complexity in customary LBM's (potentially much 
higher than that of so-called direct solution methods), which makes them inefficient or even 
potentially unsuitable for large-scale simulations in fluids. A fundamental issue is, there- 
fore, related to the construction of more accurate, or higher-order, LBM's, applicable for 
arbitrary values of the relevant physical (and asymptotic) parameters. However, the route 
which should permit to determine them is still uncertain, since the very existence of an 
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underlying exact (and non-asymptotic) discrete kinetic theory, analogous to the continuous 



inverse kinetic theory 
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271 ]. is not yet known. According to some authors 
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this should be linked to the discretization of the Boltzmann equation, or to the possible in- 
troduction of weakly compressible and thermal flow models. However, the first approach is 
not only extremely hard to implement |3l) , since it is based on the adoption of higher-order 
Gauss- Hermite quadratures (linked to the discretization of the Boltzmann equation), but its 
truncations yield at most asymptotic theories. Other approaches, which are based on 'ad 
hoc' modifications of the fluid equations (for example, introducing compressibility and/or 
„ me effects H), by deflation jno« p r ov ld e exact Navier-Sto.es solve,. WW 
critical issue is related to the numerical stability of LBM's ]8|, usually attributed to the 
violation of the condition of strict positivity (realizability condition) for the kinetic distri- 
bution function j^, 33]. Therefore, according to this viewpoint, a stability criterion should 
be achieved by imposing the existence of an H-theorem (for a review see 341]) . In an effort 
to improve the efficiency of LBM numerical implementations and to cure these instabilities, 



there has been recently a renewed interest in the LB theory. Several approaches have 



been 



proposed. The first one involves the adoption of entropic LBM's (ELBM [33] , l35l . l36l . 137 ] 
in which the equilibrium distribution satisfies also a maximum principle, defined with re- 
spect to a suitably defined entropy functional. However, usually these methods lead to 
non-polynomial equilibrium distribution functions which potentially result in higher com- 
putational complexity and lower numerical accuracy 38| . Other approaches rely on the 
adoption of multiple relaxation times. However the efficiency, of these methods is still in 
doubt. Therefore, the search for new [LB] models, overcoming these limitations, remains an 
important unsolved task. 



IV. ASYMPTOTIC LBM'S : COMPUTATIONAL COMPLEXITY 

An alternative to direct solution methods, which can reduce significantly the complexity 
caused by Poisson equation, may be achieved by so-called particle simulation methods, in 
which the dynamics of fluids is approximated in terms of a set of test particles which advance 
in time in terms of suitable evolution equations defined in such a way to satisfy identically the 
Poisson equation. Particle simulation methods rely typically on appropriate kinetic models 
for the fluid (or magnetofluid) equations which permit the evaluation of the fluid fields 
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in terms of suitable expectation values (or momenta) of the kinetic distribution function 
/(r, v, t), being respectively r and v the position an velocity of a test particle with probability 
density /(r, v, t). These kinetic models can be continuous or discrete in phase space, yielding 
respectively continuous or discrete kinetic models for the fluids. In particular, discrete 
models are those in which the kinetic distribution function is discretized in some sense 
(such as configuration space and/or velocity space and/or time), which usually leads to the 
approximate description of a continuum (the fluid) in terms of a discrete and finite set of 
test particles, each described by a suitable distribution function /j. An example is provided 
by LBM's (Lattice- Boltzmann methods) in which the distribution function /, is discretized 
in velocity space, i.e., is defined only for a suitable set of discrete velocities aj (i = 0,n) 
- to be identified with the velocities of test particles - which are all assumed as constant 
in time. However, particle simulation methods - and in particular LBM's - may exhibit, in 
their turn, different sources of undesirable algorithmic complexity. In particular, so-called 
asymptotic kinetic theories, i.e. kinetic theories which recover the prescribed set of fluid 
equations only in suitable asymptotic limits, will depend necessarily from dimensionless 
and positive parameters e±,..,ek, all assumed <C 1, and demand typically that the kinetic 
distribution function be determined accurate to prescribed order in ei, ..,£&. For example, 
the parameters e±, ..,€k my be identified, in some cases, with 1/iVfc, bein g N k the so-called 
Knudsen number, Nm the Mach number, etc.- Customary LBM's [11, 12, 13, 14, Q, [3] are 
typically constructed in such a way to satisfy the exact fluid equations (INSE) only in an 
asymptotic sense and are therefore asymptotic. These requirements inevitably give rise to 
additional computational complexity in numerical simulations based on asymptotic kinetic 
theories which are nevertheless free of the Poisson equation complexity described above. It 
is interesting to mention the possible sources of algorithmic complexity affecting customary 
LB methods, which are related to the choices of the time and space discretization scale 
lengths AL and At adopted in customary LBM's (which determine the number of cell and 
elementary time intervals in which the fluid domain and the time interval are divided). 
They are all essentially a consequence of the fact that these numerical schemes are based 
on asymptotic kinetic theories, i.e., they are characterized by an asymptotic parameter e, 
to be assumed non negative and <C 1. The parameter e enters the theory through the 
discrete kinetic distribution function f if which describes the time evolution of the i—th test 
particle in a given position (node). In particular in order that the kinetic theory recovers the 
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correct fluid equations the kinetic distribution function must remain at all times suitably 
close to an appropriate "equilibrium" distribution function f- e9 \ in the sense that denoting 
Sfi = fi — fl q the "deviation" from the equilibrium distribution, there must result 



Sfi ~ o(s) (1) 

(relaxation condition for the kinetic distribution function). In asymptotic LBM's this con- 
dition is usually satisfied by adopting an LB-BGK kinetic equation (a kinetic equation with 



a BGK collision operator [39|) which is characterized by a relaxation time r > 0. On the 
other hand, since these methods are all based on the Euler approximation for the streaming 
operator, a basic consequence [of these assumptions] is that the amplitude of the time step 
(At) used to advance in time the kinetic distribution function /j must result such that 

At 

o(e) < 1. (2) 

r 

In customary LBM's the parameter r is usually linearly proportional to the kinematic vis- 
cosity of the fluid (v). For example for the 9Q2D-(p — V)-LB scheme there results v = c 2 r/3. 
Hence it follows 

— o(e)«l. (3) 

(v— complexity of LBM's). It follows that At can become very small for weakly viscous 
fluids, with the consequence of increasing significantly the computational complexity of 
asymptotic LBM's (with respect to direct solution methods). 

Another potential source of complexity is given by the requirement that the fluid fields 
must be suitably smooth and slowly varying on the relevant discretization scales, i.e., in 
particular such that they satisfy the asymptotic ordering 

AL At 



(4) 



L 1 T 

(smallness of discretization scales), being AL and At, respectively, the size of the spatial 
cells in which the fluid is discretized and the time step for advancing in time the discretized 
kinetic distribution function, while L and T are the characteristic length and time scales 
of the fluid fields. Due to the low-order approximations used in customary LB methods 
these are expected to reproduce the correct fluid equations only in an average sense on the 
relevant characteristic scales L and T and not pointwise, i.e., in each node in which the fluid 
is discretized. This justifies, for example, the introduction of so-called "half-way bounce 
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back" boundary conditions in which the conditions of incompressibility and no-slip at a 
fixed wall in contact with the fluid are satisfied only in an average sense in the boundary 
cells. Finally, a further source of complexity for customary LB methods is that they satisfy 
the condition of incompressibility only when a suitably defined effective Mach number, M e ^ , 
results small enough, which in this case must be at least of order 

M eff ~ o(e a ) < 1, (5) 

where typically a = 3/2 (small effective Mach number assumption). In literature the pa- 
rameter M e ff is usually identified with the ratio between the magnitude of the flow velocity 
V and the constant test particle velocity c : 

V 

M eff = -, (6) 
c 

which means that the velocity of test particles must be much larger than the flow ve- 
locity, a fact which, by itself, produces a significant source of computational complex- 
ity (M e ff — complexity of LBM's). A further consequence is that the Courant number 
(Nc = ^7^, being L the magnitude of the local scale length used for the spatial discretization 
of the fluid domain) of asymptotic LBM's results of order 

N c ~ o{e p ) < 1, (7) 

with (3 > a. As a consequence the amplitude of the time step (At) used in customary 
asymptotic LBM's to advance in time the fluid fields results of order 

At ~ M e "^At 0pt (8) 

(Nq— complexity of LBM's), where Ato P t is the time step adopted by optimized numerical 
solution methods in CFD (such as those based on spectral methods), for which the Courant 
number Nc results of order unity, while AL and L denote respectively the local grid size 
adopted by the LBM and the local characteristic scale length of the fluid fields. Despite 
the absence of the N-complexity for them, asymptotic LBM's may result potentially (much) 
slower than optimized numerical methods (such as spectral methods). In fact, as indicated 
above they typically require both M e ^ 1 and ^ <C, namely At <C Ato v t- 
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V. CONCLUSIONS 



A fundamental issue in CFD is therefore the search of algorithms with reduced algorith- 
mic complexity and at the same time with improved algorithmic efficiency. This problem 
has been investigated in the framework of so-called inverse kinetic approaches, i.e., kinetic 
theories which are able to provide, with prescribed accuracy, all fluid equations expressed as 
suitable moment equations of the relevant kinetic equation, to be assumed either continuous 



40j, |4l| or discrete [42] . As indicated above, the CMFD Consortium has promoted and 



developed in the last few years an intense research activity. Part of the research activity 
has concerned inverse kinetic theories for Newtonian isothermal and non-isothermal fluids. 
In particular, first, continuous, non asymptotic, inverse kinetic theories for incompressible 
Newtonian fluids have been investigated, with the goal of developing optimal algorithms 
which do not exhibit the feature of the N— complexity and, at the same time, do not require 
subsidiary asymptotic conditions to be satisfied, such as the requirement of "closeness" in 
some sense to a suitable kinetic equilibrium. In particular, a basic features of the inverse 
kinetic theory developed is that only mild restrictions must be satisfied by the kinetic dis- 
tribution function in order to satisfy the relevant fluid equations, while the fluid pressure 
is advanced in time without solving explicitly Poisson's equation for the fluid pressure. An 
interesting issue is the possibility of formulating for incompressible Newtonian fluids a dis- 
crete kinetic approach, based on a Lattice-Boltzmann scheme which avoids or reduces the 
computational bottlenecks characteristic of previous numerical solution methods of this type 
and, at the same time, yield with prescribed accuracy the solutions of the fluid equations. 
This includes the basic feature of determining the fluid pressure without solving explicitly 
the Poisson equation. 
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